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ABSTRACT 



Aims. We present here a new method, MMF, for automatically segmenting cosmic structure into its basic components: clusters, filaments, and 
walls. Importantly, the segmentation is scale independent, so all structures are identified without prejudice as to their size or shape. The method 
is ideally suited for extracting catalogues of clusters, walls, and filaments from samples of galaxies in redshift surveys or from particles in 
cosmological N-body simulations: it makes no prior assumptions about the scale or shape of the structures. 

Methods. Our Multiscale Morphology Filter (MMF) method has been developed on the basis of visualization and feature extraction techniques 
in computer vision and medical research. The density or intensity field of the sample is smoothed over a range of scales. The smoothed signals 
are processed through a morphology response filter whose form is dictated by the particular morphological feature it seeks to extract, and 
depends on the local shape and spatial coherence of the intensity field. The morphology signal at each location is then defined to be the one 
with the maximum response across the full range of smoothing scales. 

The success of our method in identifying anisotropic features such as filaments and walls depends critically on the use of an optimally defined 
intensity field. This is accomplished by applying the DTFE reconstruction methodology to the sample particle or galaxy distribution. 
Results. We have tested our MMF Filter against a set of heuristic models of weblike patterns such as are seen in the Megaparsec cosmic matter 
distribution. To test its elfectiveness in the context of more realistic configurations we also present preliminary results from the MMF analy- 
sis of an N-body model. Comparison with alternative prescriptions for feature extraction shows that MMF is a remarkably strong structure finder 
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1. Introduction 

On scales from a few Megaparsecs up to more than a hun- 
dred Megaparsecs, the spatial cosmic matter distribution dis- 
plays a salient and pervasive weblike pattern which is per- 
ceived in the first instance as a cellular structure. The dis- 
tribution of galaxies in large scale redshift surveys such as 
the 2 -deg Field Galaxy Redshift Survey (2dF: CoHess et al 
2003 ) and the Sloan Digital Sky Survey (SPSS: IVork etal 
2000 l) clearly delineate this Cosmic Web (Ifiond et all 1 199^ 
(see van de Weygaert (2002) for a review). 

Large computer simulation s of the evolution of cosmic 
structure ("Spring el et al. , 2005 ) show prominent cellular pat- 
terns arising from gravitational instability. Galaxies accumulate 
in flattened waUs, elongated filaments and dense compact clus- 
ters. The se structures surround large near-empty void regions 
(IZeldovi ch. Einasto & Shandarin, 1982). Their spatial distribu- 



Send ojfprint requests 
miguel@astro . rug . nl 



to: 



M. 



Aragon-Calvo, e-mail: 



tion displays a distinctive frothy texture, interconnected in a 
cosmic weblike pattern. 

While it is rather straightforward to find qualitative descrip- 
tions of the spatial structure and components of the cosmic 
web, a useful, and physically meaningful, quantitative analy- 
sis has proven to be far from trivial. This would be important, 
for example, when we wish to study the effect of environment 
on the formation of galaxies and their halos. 

1.1. Multi-scale analysis 

We present here a new method for automatically segmenting 
cosmic structure into its basic components: clusters, filaments, 
and walls. Importantly, the segmentation is scale independent, 
so all structures are identified without prejudice as to their size 
or shape. 

There are two parts to this: firstly, the reconstruction of 
a continuous density field from a point sample and secondly, 
the identification of structures within that density field. For 
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the first part we use t he Delaunay Tessellation Fie l d Esti mator 
(DTFE) technique of ISchaap & van de Wevgaerl j2000l) . The 
second part, which is the main thrust of this paper, consists of 
a series of morphology filters that identify, in a scale indepen- 
dent manner, particular kinds of structure in data. The method 
is referred to as the Multiscale Morphology Filter (MMF) and 
is based on the kind of Scale Space analysis that has in recent 
years proved so successful in imaging science. 

It is worth emphasising at this juncture that we have chosen 
a specific implementation of this kind of multi-scale analysis. 
Our choice is made on the following grounds: (a) it is simple 
to understand and program, (b) it works under quite general 
conditions and (c) the approach is generic and easy to mod- 
ify. There are many alternative multi-scale strategies: we leave 
those for another day or for other people to follow up. Thus we 
shall try to keep this presentation as general as possible so that 
the points at which we make implementation specific choices 
are clear 



2. Structure finding 

Many attempts to describe, let alone identify, the features and 
components of the Cosmic Web have been of a mainly heuris- 
tic nature. There is a variety of statistical measures character- 
izing specific aspects o f the large scale matter d istribution (for 
an extensive review see lMartmez & Saai , l2002h . For complete- 
ness and comparison, we list briefly a selection of methods for 
structure characterisation and finding. It is perhaps interesting 
to note two things about this list: 

a) each of the methods tends to be specific to one particular 
structural entity 

b) there are no explicit wall-finders. 

This emphasises an important aspect of our Scale Space ap- 
proach: it provides a uniform approach to finding Blobs, 
Filaments and Walls as individual objects that can be cata- 
logued and studied. 



1.2. Emergence of hierarchical web-like structure 



2.1. Structure from higher moments 



Structure in the Universe emerged as a result of the gravita- 
tional growth of small amplitude primordial density and ve- 
locity perturbations. Following the initial linear growth of the 
Gaussian primordial perturbations, the gravitational clustering 
process leads to the emergence of complex patterns and struc- 
tures in the density field. At least three characteristics of the 
midly nonlinear cosmic matter stand out. 

The most prominent property is its hierarchical nature. The 
gravitational clustering process proceeds such that small struc- 
tures are the first to materialize and subsequently merge into 
ever larger entities. As a result each emerging cosmic structure 
consists of various levels of substructure. Hence, upon seeking 
to identify structure at one characteristic spatial scale we need 
to take into account a range of scales. 

The second prominent aspect is that of the weblike geome- 
try marked by highly elongated filamentary and flattened planar 
structures. The existence of the cosmic web can be understood 
through the tendency of matter concentrations to contract and 
collapse gravitationally in an anisotropic manner. 

A final conspicuous aspect is that of the dominant presence 
of large roundish underdense regions, the voids. They form in 
and around density troughs in the primordial density field. 

The challenge for any viable analysis tool is to trace, high- 
light and measure these features of the cosmic web. 

1.3. Outline of this paper 

We start in section [3] by reviewing the DTFE method that is 
used to sample discrete point sets onto a regular mesh. Then 
in section|5]we introduce the basic ideas from scale space the- 
ory that we will use. In section |5] we introduce the morphology 
filters and give them a geometrical interpretation. The filters 
are tested using a Voronoi model in section|9] We present brief 
results from an N-body simulation in section [10] leaving a de- 
tailed study to a subsequent paper in this series. 



The clustering of galaxies and matter is most commonly 
described in terms of a hierarchy of correlation functions. 
The two-point correlation function (and its Fourier trans- 
form, the power spectrum) remains the mainstay of cos- 
mological clustering analysis and has a solid physical ba- 
sis. However, the nontrivial and nonlinear patterns of the 
cosmic web are mostly a result of the phase correlations 



in the cosmic matter distribution ([Rvden & Gramannl II991 
IChiang & Coiesil200(]HColes & Chiangll2000l) . Whil e fliis i n 



formation is contained in the moments of cell counts ("Peebles', 
1980; de Lapparent, Gefler & Huchra, 1991; Gaztanaga, 1992) 
and, more formally so, in the full hierarchy of M-point corre- 
lation functions ^m, their measurement has proven t o be im- 
practical for all but the l owest orders (Peebles. 1980;,Szapudi , 
ll998llJones etallbOOSh . 

The Void pr obability Function ('Whitel Il979 : 
iLachieze-Rey. da Costa & Mauro gordato.. 1992 ) provided 



a characterisation the "voidness" of the Universe in terms 
of a function that combined information from many higher 
moments of the point distribution. But, again, this has not 
provided any identification of individual voids. 

2.2. Topological methods 

The shape of the local matter distribution may be traced on the 
basis of an analysis of the stat istical properties of its i nertial 
moments ( Babul & Starkmanl 1992 : Luo & Vishniad 1 1995 



Basilakos. Plionis & Rowan-Robinsoni 120011) . These concepts 



are closely related to the full characterization of the topology 
of the matter distribution in terms of four Minkowski func- 



tional s (iMecke. Buchert & Wagnerl Il994t ISchmalzing et al 
1999b . They are solidly based on the theory of spatial statis- 
tics and also have the great advantage of being known an- 
alytically in the case of Gaussian random fields. In partic- 
ular, the genus of the density field has received substan- 
tial attention as a strongly discriminating factor between in- 
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trinsicaUy different spatial patt erns dGott. Dickinson & Melott , 
1986HHovle & VogelevlEool . 

The Minkowski functionals provide global charac- 
terisations of structure. An attempt to extend its scope 
towards providing locally defined topological mea- 
sures of the density field has been developed in the 
SURFGEN projec t defined by Sahni and Shandarin and 



maxima and saddle points. This is computationally intensive 
but quite eff'ective, though it does depend on the artefacts in 
the reconstruction of the continuous density field. 

Stoica et al. (l2005h use a generalization of the classical 



their coworkers dSahni Sathyaprakash & Shandarin , 1998 : 
Shandarin. Sheth & Sahnil 120041) . The main problem remains 
the user-defined, and thus potentially biased, nature of the 
continuous density field inferred from the sample of discrete 
objects. The usual filtering techniques suppress substructure 
on a scale smaller than the filter radius, introduce artificial 
topological features in sparsely sampled regions and diminish 
the flattened or elongated morphology of the spatial patterns. 
Quite possibly the introduction of more advanced geometry 

based methods to trace the density field may prove a major ^002; Plionis & Basilakos, 
advance towards solving this problem. 



Candy model to locate and catalogue filaments in galaxy sur- 
veys. This approach has the advantage that it works directly 
with the original point process and does not require the creation 
of a continuous density field. However, it is very computation- 
ally intensive. 



2.5. Void Finding 

Voids are distinctive and striking features of the cosmic 
web, yet finding them systematically in surveys and simu- 
lations has proved rather difficult. There have been exten- 
sive searches for voids in g alaxy catalogues (IHoyle & Vogeley . 



jArbabi-B idgoh & Miillei 



2002 ) and in numerical simulations 



2002HAikio & Mahonen[ll998h . 



Importantly, [Martinez et al.1 ( l2005h and ISaar et all (120070 
have generalized the use of Minkowski Functionals by cal- 
culating their values in a hierarchy of scales generated from 
wavelet-smoothed volume limited subsamples of the 2dF cata- 
logue. This approach is particularly effective in dealing with 
non-Gaussian point distributions since the smoothing is not 
predicated on the use of Gaussian smoothing kernels. 

2.3. Cluster finding 

In the context of analysing distributions of galaxies we can 
think of cluster finding algorithms. There we might define a 
cluster as an aggregate of neighbouring galaxies sharing some 
localised part of velocity space. Algorithms like HOP attempt 
to do this. However, there are always issues arising such as how 
to deal with substructure: that perhaps comes down to the defin- 
tion of what a cluster is. Here we focus on defining coherent 
structures based on particle positions alone. The velocity space 
data is not used since there is no prior prejudice as to what the 
velocity space should look like. 

2.4. Filament finding 

The connectedness of elongated supercluster structures in the 
cosmic matter distribution was first probed by means of perco- 
lation anal ysis, introduced and emphasized by Zel 'dovich and 
coworkers dZeldovich. Einasto & Shandarinill982h . while a re- 
lated graph-theoretical construct, the minimum spanning tree 
of the galaxy distribution, was ext ensively probed and analysed 
by Bhavsar and collaborators dBarrow. Bhavsar & Sonoda , 
1985l;lGrah"amlll995l;[Colberd,l2007b in an attempt to develop 
an objective measure of filamentarity. 

Finding filaments joining neighbouring clusters 
h as been tackled, using qu it e difi 'erent techni ques, by 



Pimbblet 



Colberg. Krughofi^ & ConnoUvl ( 120051) and by . 
( 2005b . More general filament finders have been put forward 
by a number of authors. Skeleton an alysis of the density 
field dNovikov. Colombi & Pore , l2006l) describes continuous 
density fields by relating density field gradients to density 



Several factors contribute to making systematic void- 
finding difficult. The fact that voids are almost empty of galax 
ies means that the sampling density p lays a key rol e in de 



termining what is or is not a void dSchmidt et all 1200 Ih . 



Moreover, void finders are of ten predicated on bui l ding v oid 
structures out of cubic cells ( KaufFma nn & FairallL 1991) or 
out of spheres (e.g: Patiri et al.- .2006.) . Such methods attempt 
to synthesize voids from the intersection of cubic or spher- 
ical elements and do so with varying degrees of success. 
The Aspen-Amsterdam Void Finder Comparison Project of 



Colberg. Pearce et al.l d2007h will clarify many of these issues. 

The Watershed-based algorithm of 

Platen, van de Wevgaert & JonesI d2007b aims to avoid is- 



sues of both sampHng density and shape. 
2.6. Structure from Scale Space 

Combining the local Hessian matrix eigenvalues on various 
scales, this is the new technique that we are presenting here 
for the first time in the cosmological context. 

Scale space analysis looks for structures of a mathemati- 
cally specified type in a hierarchical, scale independent, man- 
ner It is presumed that the specific structural characteristic is 
quantified by some appropriate parameter (e.g.: density, eccen- 
tricity, direction, curvature components). The data is filtered to 
produce a hierarchy of maps having different resolutions, and 
at each point, the dominant parameter value is selected from the 
hierarchy to construct the scale independent map. We refer to 
this scale-filtering processes as a Multiscale morphology filter . 

For simplicity, the paper describes one specific implemen- 
tation, or embodiment, of the process in relation to the problem 
of cataloguing the structural elements of the cosmic web. Other 
embodiments are possible, but the present one turns out to be 
highly effective in structure segregation and feature identifica- 
tion. 

While this sounds relatively straightforward, in practise a 
number of things are required to execute the process. Firstly 
there must be an unambiguous defintion of the structure- 
defining characteristic. In the present case we shall use the prin- 
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cipal components of the local curvature of the density field at 
each point as a morphology type indicator. This requires that 
the density be defined at all points of a grid, and so there must 
be a method for going from a discrete point set to a grid sam- 
pled continuous density field. We choose to do this using the 
DTFE methodology since that does minimal damage to the 
structural morphology of the density field. 

Since we are looking for three distinct structural morpholo- 
gies, blobs, walls and filaments, we have to apply the segmen- 
tation process three times. However, since we shall be using 
curvature components as structural indicators, we shall have to 
eliminate the blobs before looking for filaments, and we shall 
then have to eliminate the filaments before looking for walls. 



3. Resampling and Rescaling Point sets 

The cosmological problem presents its own difficulties, not the 
least of which is the fact that the data set is presented, not as 
a density field, but as a set of discrete points which are pre- 
sumed to sample some underlying density field. However, the 
filtering procedures we use here for defining objects act on con- 
tinuous fields (or images) and require continuous first and sec- 
ond derivatives of field values. It is therefore necessary to re- 
sample the point set data on a grid. In doing this we need to 
assure ourselves that the objects, structures, features and pat- 
terns in these fields are resampled in an optimal way: both sub- 
structure and morphological characteristics must be preserved. 
To achieve this we use the DTFE reconstruction of the density 
field. 



3.1. The DTFE density field 

The Delaunay Triangulation Field Estimato r ("DTFE") 



dSchaap & van de Wevgaert , l2000t l ISchaapll2007h is a powerful 
new method, based upon concepts from computational geom- 
etry ( Okabe et al. . 2000l) that ofifers a "safe" reconstruction in 
that it accurately preserves the local features. DTFE produces 
a morphologically unbiased and optimized continuous density 
field retaining all features visible in a discrete galaxy or particle 
distribution. 

The input samples for our analysis are mostly samples 
of galaxy positions obtained by galaxy redshift surveys or 
the positions of a large number of particles produced by 
N-body simulations of cosmic structure formation. In or- 
der to define a proper continuous field from discrete dis- 
tribution of points - computer particles or galaxies - we 
translate the spatial point sample into a continuous density 
field by means of the Delaunay Tessellation Fie ld Estimator 
Schaap & van de Weygaerl (l200d) : [Scha^ (l2007h . 

The DTFE technique recovers fully volume-covering and 
volume-weighted continuous fields from a discrete set of 
sample field values. The method has be e n dev eloped by 



(iSchaap & van de Weygaerl . 20001: I Schaapl l2007l ,see) and 



forms an el aboration of the velocity interpolation scheme in- 
troduced by lBernardeau & van deWevgaertl ( ll996l) . It is based 
on the use of the Voronoi and Delaunay tessellations of a given 
spatial point distribution. It provides a basis for a natural, fully 



self-adaptive filter in which the Delaunay tessellations are used 
as multidimensional interpolation intervals. 

The primary ingredient of the DTFE method is the 
Delaunay tessellation of the particle distribution. The Delaunay 
tessellation of a point set is the uniquely defined and volume- 
covering tessellation of mutually disjunct Delaunay tetrahedra. 
A Delaunay tetrahedron is defined by the set of four points 
whose circumscribing sphe re does not contain any of the other 
points in the generating set ('Delaunay , 1934 ) (triangles in 2D). 
The Delaunay tessellation is intimately related to the Voronoi 
tessellation of the point set, they are each others dual. The 
Voronoi tessellation of a point set is the division of space into 
mutually disjunct polyhedra, each polyhedron consisting of the 
part of space closer to the defining point than any of the other 



pomts (YoTonm, 



1908: Okabe I 



inmg pomt tr 
et al.L l2000^. 



DTFE exploits three particular properties of Voronoi and 
Delaunay tessellations. The tessellations are very sensitive to 
the local point density. The DTFE method uses this fact to de- 
fine a local estimate of the density on the basis of the inverse of 
the volume of the tessellation cells. Equally important is their 
sensitivity to the local geometry of the point distribution, which 
allows them to trace anisotropic features such as encountered in 
the cosmic web. Finally it uses the adaptive and minimum tri- 
angulation properties of Delaunay tessellations to use them as 
adaptive spatial interpolation intervals for irregular point dis- 
tributions. In this it is the first order version of the Natural 

I 1 r 



Nei}>hbour method (NN method: Sibson, 1980, 1981; Wats^, 
1992; Braun & Sambridge..l995:.Sukumar. 1998;.Okabe etal 



2000). 

One of the important - and crucial - properties of a pro- 
cessed DTFE density field is that it is capable of delineating 
three fundamental characteristics of the spatial structure of the 
Megaparsec cosmic matter distribution. It outlines the full hi- 
erarchy of substructures present in the sampling point distribu- 
tion, relating to the standard view of structure in the Universe 
having arisen through the gradual hierarchical buildup of mat- 
ter concentrations. DTFE also reproduces any anisotropic pat- 
terns in the density distribution without diluting their intrinsic 
geometrical properties. This is a great advantage when seek- 
ing to analyze the cosmic matter distribution, characterized 
by prominent filamentary and wall-like components linking up 
into a cosmic web. A third important aspect of DTFE is that 
it outlines the presence and shape of voidlike regions. DTFE 
renders the low-density regions as regions of slowly varying, 
moderately low density values through the interpolation defini- 
tion of the DTFE field reconstruction. 

An outline of the DTFE reconstruction procedure can be 
found in appendixlAl 

3.2. Rescaling 

In building the scale space we need to construct a hierarchy 
of rescaled replicas of the original grid-sampled data. In this 
paper this is done simply by applying a hierarchy of isotropic 
Derivative of Gaussian smoothing filters to the data. 

Of course, substructure and morphological characteristics 
will be altered during this hierarchical smoothing process. The 
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Fig. 1. DTFE image of a slice through the N-body simulation used in this work. Left: DTFE density field in a central slice. Right; 
the corresponding particle distribution in a slice of width 5/? 'Mpc. 



smearing of features through smoothing is inevitable if we 
smooth using isotropic filters and there has been some dis- 
cussion as to whether one might do better by rescaling in 
such a way as to minimise feature smearing (for example 



Martinez et a l.. 2005; Saar et al., 2007). It is possible to use re- 



fined (nonlinear) smoothing procedures that minimize the side 
effects of smoothing but that issue is not addressed here. Here, 
we simply rescale using isotropic Gaussian filters; this seems 
to work very well and avoids complications arising from using 
other filters. 



4. Scale Space Analysis 

In this contribution we introduce a method for recognizing and 
identifying features in data b ased on the use of a "Scale Space" 



repr esentation of the data ( Florack et al. . 1992 : Lindeber^ 



repr e 
[1998 



I199&) . The Scale Space representation of a data set consists 
simply of a sequence of copies of the data having different res- 
olutions. A feature searching algorithm is applied to all of these 
copies, and the features are extracted in a scale independent 
manner by suitably combining the information from all copies. 

We use a particular feature recognition process based 
on eigenvalues of the Hessian matrix of the density field. 
It should be understood that the technique we describe 
here could well be used with other feature recognition 
systems, such as, for example, t he Sh apeFinder process 



dSahni. Sathyaprakash & ShandarinL Il998l) . Scale Space is a 
powerful tool for scale independent data analysis. 

4.1. Image processing 

The use of this technique can be tra ced back to the work o f 
David Marr at M.I.T in the 1970's dMarr & Hildrethi Il998h . 



reviewed in his seminal book on the physiology of image un- 
derstanding; Vision (Marr, 1980). There (loc. cit. Chapter 2, 
especially figures 2-10 and 2-23) he describes what is called 
the "Primal Sketch" and the use of what today are called "Man- 
Wavelets" in extracting scale independent information. We ap- 
ply precisely this transformation to a scale space representation 
of a cosmological density field, and in doing so ostensibly ex- 
tract features in much the same way, according to Marr, that the 
human visual cortex does. 



More recently, Fran gi et alJ (119981) and ISato et all (119981) 
used Scale Space analysis for detecting the web of blood ves- 
sels in a medical image. The vascular system is a notoriously 
complex pattern of elongated tenuous features whose branch- 
ing make it closely resemble a fractal network. We translate, 
extend and optimize this technology towards the recognition of 
the major characteristic structural elements in the Megaparsec 
matter distribution. The resulting methodology yields a unique 
framework for the combined identification of dense, compact 
bloblike clusters, of the salient and moderately dense elongated 
filaments and of tenuous planar walls. 

4.2. Multiscale Structure Identification 

Segmentation of a density field into distinct, meaningful, com- 
ponents has been one of the major goals of image processing 
over the past decades. There are two stages involved; firstly 
providing a criterion describing the basis for the segmenta- 
tion, be it colour, texture, motion or some other attribute and 
secondly providing an algorithm whereby those distinguish- 
ing attributes can be automatically and unambiguously iden- 
tified. Ambiguities in structure finding frequently occur when 
the sought-for structure exists on a variety of scales that may 
be nested hierarchically. 
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Symbol 


Name 


Description 


Eqn 




Scale Space Map 


Combination filtered density maps /s,n over all levels n. 


(0 


£ 
S 


Morphology Mask 
Shape Significance Map 


Region of space obeying shape constraint. 

F — 1 ■ Inpfitinir; nhpvinQ 'jhanp pnTTjtvaint 

E=0: locations not obeying shape constraints 
Feature shape fidehty for each point locale. 

A/fpaQiirp'; rnnfnnnfinr'P tn Inral ^hanp mtpria 


(l6l 
(l9l 


M 


Morphology Map 


Soft thresholded version of <S. The threshold selects out the most 

locally shape conformant features. Requires input of a threshold parameter j8 


Go) 


I 


Morphology Intensity Map 


Map of /^3 for blobs, A2 for filaments or Ai for walls 

Modulates Morphology map, meant to avoid enhancing noisy low intensity structures 




T 
1 


iviui uinjiug y ni LCI 


V^UilSLl U.L.LCQ ilUIIl-f dilU yrl. IVlUi UllUlUgV WClglllCU. IIILCI 

for the Morphology Mask. Provides each location which obeys the morphology constraint 
with a measure of the strength of morphology signal. 


ca 


T 


Feature Map 


Product of morphology mask £ and corresponding morphology filter T. 

There is one Feature Map for each level in the Scale-Space, representing local structures as 

seen on the different scales of the Scale-Space 






Scale-Space Map Stack 


Constructed from the fj for all levels in the Scale-Space 

Each pixel in this map is the greatest value of the corresponding pixels 

in the Feature maps that make up the Scale-Space stack 


Gl 





Object Map 


Inclusion of astrophysical & cosmological criteria to select physically recognizable objects 
Produced by thresholding Scale-Space Map Stack 

Threshold criterion determined by cosmological/astrophysical considerations 


(sec.O 



4.3. the Multiscale Morphology Filter: Outline 



Each pass involves the following components and procedures: 



The technique presented here, the Multiscale Morphology 
Filter (MMF), looks to synthesize global structures by identify- 
ing local structures on a variety of scales and assembling them 
into a single scale independent structural map. The assembly is 
done by looking at each point and asking which of the struc- 
tures found on the various search scales dominates the local 
environment. This is the essence of the so-called Scale Space 
approach. We first provide an outline of the various stages in- 
volved with the MMF method. In the subsequent sections we 
treat various aspects in more detail. 

4.4. The Analysis Cycle 

We are looking for three distinct morphologies within the same 
distribution. This requires three passes through the data, each 
time eliminating the features found in the previous pass. In 
the first pass, the blobs in the dataset are identified along with 
their enclosed datapoints. The points that are in blobs are elimi- 
nated and then the filaments are identified with their constituent 
points. After eliminating the filament points the walls and their 
constituent points can be identified. 



Point Dataset 

For each pass this is the set of galaxies or particles in 
an N-body model from which we are going to extract a 
specified feature. In the first pass this is the full data sample 
within which we are going to identify blobs. On the sec- 
ond pass it is the original point set from which the points 
in the blobs have been removed. Likewise for the third pass. 

DTFE Density Field 

The discrete spatial distribution of galaxies, or particles in a 
N-body computer model, is resampled to give a continuous 
volume-filling density field map /dtfe on a high resolution 
grid. In order to guarantee an optimal representation of 
morphological features this is accomplished on th e basis 
of the DTFE method dSchaap & van de Wevgaerl I2OOOI; 
Schaapil2007l) . 



Scale filtering 

The raw DTFE density field /dtfe is filtered over a range 
of spatial scales in order to produce a family O of 
smoothed density maps fs , n, each defining a level of the 
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Scale-Space representation. The range of scales is set by 
the particular interest of the analysis. 

Hessian & Eigenvalues 

The Hessian matrix V^/s of the density field is computed 
at each point of each of the smoothed density fields in 
the filtered Scale-Space density maps fs- At each point 
the eigenvalues Ak (k - 1,2,3) of the Hessian matrix are 
determined. 

Morphology Mask 

The Morphology Mask fimorph identifies the locations 
obeying the required morphology/shape constraints. At 
every location in every map, £ = 1 if the shape constraint 
is valid, fi = if it does not. This is a "hard" filter. 

Shape Significance Map 

A Feature shape Significance (or fidelity) index »Smorph is 
determined for the specified morphology. This is done on 
the basis of the signs and ratios of the three eigenvalues 
Ak {k = 1,2,3), and thus dependent only on the local 
variations of the field on the various scales present in the 
scale space maps 

Morphology Response Map 

The Morphology Response Filter, Almm-phi is the soft 
thresholded version of the shape significance map >Smorph- 
It selects out the most locally shape conformant features 
and is computed for each scale space level by processing 
>Smorph, weighted by a specified threshold parameter /3. 

Morphology Intensity Map 

In order to avoid enhancing noisy low intensity structures 
we include a Morphology Intensity function Xmoiph that 
modulates the morphology response map according to 
some measure of the feature strength. We characterise 
feature strength by the values of the specifc eigenvalues: 
Al for the walls, A2 for the filaments and A3 for the blobs. 

Morphology Filter 

Morphology weighted filter T^noiph for the Morphology 
Mask £morph- Provides each location which obeys the 
morphology constraint with a measure of the strength of 
morphology signal. 

Feature Map 

For each level of Scale-Space the feature map 'Fmorph 
is constructed from the Feature Intensity Map Jmorph 
and the Morphology Response Map. This represents local 
structures as seen on the different scales of the Scale-Space. 

Scale-Space Map Stack 

By combining the individual Feature Maps TL,morp\\ of 
each level of Scale-Space, the ultimate scale independent 
map of features is produced, the Scale-Space Map Stack 
Each pixel in this map is the maximum value of the 
corresponding pixels in the Feature maps that make up the 



Scale-Space stack. 

• Object Maps 

Astrophysical and Cosmological criteria determine the 
final Object Maps <9morph- These maps are produced by 
thresholding the Scale-Space Map Stack 'I' according 
to a criterion translating a feature map of physically 
recognizable objects. 

• Datapoint identification 

Datapoints within the feature contours of the object map 
<5moiph are identified. They are removed from the orig- 
inal dataset at each pass through the feature finding process. 



5. Scale Space Technology 

5.1. Scale-Space Filtering 

The so-called Scale-Space approach to morphology consists 
simply of calculating and comparing morphology indicators on 
a variety of scales. Fundamental in this is the ability to view 
a given dataset on different scales. This task is accomplished 
simply by convolving the original data f(x) with smoothing 
filters W to produce a smoothed field fsix): 



fs(x) 



dyf(y) W(y,x) 



The smoothing filter could be any of a number of suitable fil- 
ters: it is usual, though neither necessary nor optimal, to choose 
filters based on Gaussian functions. There are alternatives to 
this scaling strategy: any form or pyramidal or wavelet trans- 
form will have a similar effect. 

In this paper we generate scaled representations of the data 
by repeatedly smoothing the DTFE reconstructed density field 
/dtfe with a hierarchy of spherically symmetric Gaussian fil- 
ters Wc having different widths R: 

fs(x) ^ J /dtfeCv) WG(y,x) 

where Wc denotes a Gaussian filter of width R: 

\y-x' 



WG{y,x) 



1 

(2;r/?2)3/2 



exp 



2/?2 



(1) 



A pass of the smoothing filter attenuates structure on scales 
smaller than the filter width. 

The scale-space MMF analysis described in this study 
involves a discrete number of -H 1 levels, n - Q, . . . ,N . 
Following dSato et al.L 119981) we use a nested hierarchy of fil- 
ters having widths differing by a factor of V2: 



R„^{^)"Ro 



(2) 



The base-scale Rq is tak en to be equal to the pixel scale of the 
raw DTFE density map. Sato et al.l (11998.) showed that using a 
ratio of ~ V2 between discrete levels involves a deviation of a 
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mere 4% with respect to the ideal case of a continuum of scale- 
space levels. Q As a retrospective on this research we would 
argue that, in the context of cosmic structure, the factor of V2 
is somewhat too coarse. 

The largest structure that survives this process is deter- 
mined by the effective width of the filter used in the final 
smoothing stage. For our purposes it is sufficient to use n - 5. 

We shall denote the n''' level smoothed version of the DTFE 
reconstructed field /dtfe by the symbol f„ . 

The Scale Space itself is constructed by stacking these var- 
iously smoothed data sets, yielding the family O of smoothed 
density maps /„ : 



0) 



(3) 



levels n 



A data point can be viewed at any of the scales where scaled 
data has been generated. The crux of the concept is that the 
neighbourhood of a given point will look different at each scale. 
There are potentially many ways of making a comparison of the 
scale dependence of local environment. We chose here to use 
the Hessian Matrix of the local density distribution in each of 
the smoothed rephcas of the original data. 

5.2. The Hessian 

At each point of each dataset in the Scale Space view of the 
data we can quantify the local "shape" of the density field in 
the neighbourhood of that point by calculating, at each point 
the eigenvalues of the Hessian Matrix of the data values. 

We can express the local variations around a point xq of the 
density field f(x) as a Taylor expansion: 



where 

/aa fyx fz 

= U fyy U (5) 

^ fxz fyz fzz , 

is the Hessian matrix. Subscripts here denote partial derivatives 
of / with respect to the named variable. There are many possi- 
ble algorithms for evaluating these derivatives. 

In our case we compute the scale-space Hessian matrices 
for each level n directly from the DTFE density field, via the 
convolution 



52 



dxidxj 



fs{x) - /dtfe 



dxjdxj 



{xi - yi){xj - y,) - 5iiR\ 



(6) 



where xi,X2,X2 - x, y, z and Sjj is the Kronecker delta. In other 
words, the scale space representation of the Hessian matrix for 
each level n is evaluated by means of a convolution with the 
second derivatives of the Gaussian filter, also known as the 
Marr (or, less appropriately, "Mexican Hat") Wavelet. 

In order to properly compare the values of the Hessian aris- 
ing from the differenlty scaled variants of the data that make up 
the Scale Space we must use a renormalised Hessian: 

(7) 

where Rs is the filter width that has been used ( V2"/?o for level 
n in our case). Instead of using this 'natural' renormalization, 
it would be possible to use a scaling factor 7?^''. Using values 
y > 1 will give a bias towards finding larger structures, while 
values y < 1 will give a bias towards finding smaller structures. 



f(xo + *) = f(xo) + s^Vf(xo) + ls^mxo)s + ... (4) 



' It is interesting to note also that iMarj ( Il98(]h had already com- 
mented on the importance of the V2 factor on psycho-visual grounds. 



5.3. Eigenvalue and Eigenvectors 

The eigenvalues of the Hessian matrix evaluated at a point 
quantify the rate of change of the field gradient in various di- 
rections about each point. The eigenvalues are coordinate inde- 
pendent measures by the components of the second derivatives 
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Fig. 3. Maps of the 6 independent components of the (filtered) density field Hessian, ordered by their identity. Top row: Vn/; 
Central row: V21/, V22/; Bottom row: V31/, V32/, V33/ 



of the field at each point xq. A small eigenvalue indicates a low 
rate of change of the field values in the corresponding eigen- 
direction, and vice versa. 

We denote these eigenvalues by Aaix) and arrange them so 
that/li > /I2 > A3: 



dxidxj 



■ Aa(x) 6ij 

with A\ 



0, 



> A2> Ai 



1,2,3 



(8) 



The Aiix) are coordinate independent descriptors of the be- 
haviour of the density field in the locality of the point x and 



can be combined to create a variety of morphological indica- 
tors. The corresponding eigenvectors show the local orienta- 
tion of the morphology characteristics. Note, however, that in 
this study we do not make use of the eigenvectors. 

6. Scale-Space Feature Detection and Extraction 

The eigenvalues of the Hessian therefore encode the local mor- 
phology of the density field in terms of the curvature compo- 
nents of the local density field in the direction of the corre- 
sponding eigenvectors. Evaluating the eigenvalues and eigen- 
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Fig. 4. Maps of the eigenvalues of the Hessian matrix at 3 different scales (levels). From top to bottom: the 3 eigenvalues Ai, A2 
and A3 {Al > A2 > A3). From left to right: 3 different scales Ri Rs and/?5, {Ri > R2> R5). Positive values are represented as gray 
shades in logarithmic scale while negative values are indicated by contour lines also in logarithmic scale. 



vectors for the renormalised Hessian of each dataset in a 
Scale Space shows how the local morphology changes with 
scale. 

With the local curvature and shape encapsulated in the three 
eigenvalues Ai, A2 and Aj of the Hessian, the MMF seeks to 
identify the regions in space which correspond to a certain mor- 
phology and at the scale at which the corresponding morphol- 
ogy signal attains its optimal value. First we set out to select 
these regions by means of a Morphology Mask. Subsequently 
we develop a filter-based procedure for assigning at each scale 



a local weight which is used to select the scale at which the 
morphology reaches its strongest signal. 

6. 1 . Morphology Mask: B 

Locally "spherical" topology is indicated by all three eigenval- 
ues being similar in size, and locally "filamentary" topology 
is indicated by having two similar eigenvalues and a negligi- 
ble third; the direction of the filamentary structure is then in 
the direction of the eingenvector corresponding to the small- 
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Fig. 5. Morphology Mask £: on the basis of the 3 eigenvalues Ai , A2 and /I3 at each location we determine whether the morpho- 
logical criterion - here whether it corresponds to a filament (table 1) - is valid. If so £ = 1, otherwise it is £ = 0. Top row: maps 
of the three eigenvalues; bottom row: the Morphology Mask S. 



Structure 


A ratios 


A constraints 


Blob 


A 


^A2^Ai 


A3 


< 


/Ij < ; Al <0 


Line 


^1 


/i2 » /Is 


A3 


< 


A2<0 


Sheet 


Ai 


» /I2 ^ /Is 


A3 


< 





Table 1. Eigenvalue relationships defining the characteristic 
morphologies. The /l-conditions describe objects with inten- 
sity higher than their local background as clusters, filaments 
or walls. For voids we would have to reverse the sign of the 
eigenvalues. 



est (insignificant) eigenvalue. A locally "sheet-like' structure 
is characterised by one dominant eigenvalue, its correspond- 
ing eigenvector indicating the normal to the sheet. The formal 
morphology conditions are listed in table |6] 



There are many ways of using the eigenvalues of the 
Hessians in the Scale Space representation of the data to iden- 
tify and demarcate specific types of structure. Here we start by 
defining a morphology mask. The Morphology Mask £morph is 
a hard filter which identifies all pixels obeying the morphology 
and shape condition: 



^morph ~ 



1 morphology constraint valid 

morphology constraint invalid 
See figure|5]to see how this works. 

6.2. Feature shape fidelity: S 

The degree of "blobbiness", "filamentariness" or "wallness" is 
reflected in the degree to which the inequalities of table|6]defin- 
ing those structures are satisfied. We would be impressed by a 
blob in which all three eignevalues were equal - it would look 
like a spherical lump. We would be less impressed if there was 
a factor 3 between the eigenvalues since the blob would then 
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Fig. 6. Morphology Filter T. The Morphology Response function M (top centre) is the soft thresholded version of the Shape 
Significance map S (left frame), determined from the values of the eigenvalues Ai, A2 and /I3. The Morphology Intensity function 
I (bottom centre) is also computed from the /I's using equation ( fTTT ). Finally, the Morphology Filter T (right frame) is obtained 
by combining M with I. 



look more like a flattened sausage while not manifestly being a 
filament or a wall. 

The following shape indices reflect the strength S of the 
classification in terms of the local geometry as characterised 
by the A's. 



^morph 



Blob 



Filament 



,('-m)('-ra) 



(9) 



It is important to emphasise when using this equation that the 
values of S are only meaningful if the relevant inequalities in 
table|6]are already satisfied. 

As a cautionary warning it must be stressed that we cannot 
identify a point as being part of a locally filamentary structure 
and assess the significance by using an evaluation of S that ap- 
plies to blobs or walls. Likewise the value of S cannot be used 
to assess the relative significance of different types of struc- 
ture. This means that the identification of structural elements 
using this eigenvalue classification scheme must be done cycli- 
cally; first find blobs (three roughly equal eigenvalues), then 
lines (two roughly equal and dominant eigenvalues) and finally 



walls (one dominant eigenvalue). There are other schemes that 
are one-pass classifiers. 

We shall use the symbols >Sbiob, >Sfiiament, >Swaii to denote the 
values of S computed for each kind of feature. 

6.3. Morphology Response Filter: M 

We shall need a filter that preferentially selects out points where 
the value of the feature shape parameter S lies above some 
threshold. With this we can tune the aggressiveness of feature- 
selection. This can be done by defining a morphology measure 
Mhy 

M.»ph = l-exp(-:^) (10) 

where morph - (blob, filament, or wall). The adjustable pa- 
rameter y6 tunes the discrimination level of the morphology re- 
sponse filter. A typical value is y6 = 0.5. Lower values will 
increase the feature selectivity. Higher values will decrease the 
selectivity giving feature images with smooth features but con- 
tamination from other morphologies. 

We shall use the symbols Albiob, Alfiiament, A^waii to denote 
the values of M computed for each kind of feature. 

Methods of thresholding image data such as equation ( fTOl l 
are generally referred to as "soft thresholding", as opposed to 
"hard thresholding" in which all values below a critical value 
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Fig. 7. The Feature Map T (righthand frame) is computed for each scale and is equal to the Morphology Filter T at the locations 
where the Morphology Mask & is unity (and nonzero). 



are zeroed. Soft thresholding results in visually more appealing 
density distributions. See figure ©. 



6.4. Morphology Intensity Map I 

Morphology Intensity is a property of structures that represents 
how strong the feature is: a filament that is nice and narrow 
is in some sense more filament-like than one which is rather 
wide and diffuse. The discriminating factor in this case is the 
magnitude of the eigenvalue /I2. Note that it would be inappro- 
priate to normalise or non-dimensionalize this relative to some 
local values such as the sum of the local A's: it is the fact of 
comparing the A values at different spatial locations that dis- 
criminates features. If, in our example, the value of A2 were 
roughly constant over the data set, we would not be impressed 
by any filamentariness. 

Oian. Sone & Doi (l2003h noted that the smallest eigenvalue 



(A3) will be large only for blobs, while A2 will be large for 
blobs and filaments, and Ai for blobs, filaments, and walls. 
Combining these relations with the A constraints in table |6] we 
can use the following intensity function: 





■^3 


Blob 




^ moiph ~ 




Filament 


(11) 






WaU 





The use of this morphology intensity function solves the prob- 
lem of detecting low-intensity/noisy structures but it introduces 
another problem: the range of values of Imorph is not well de- 
fined within a given interval since it depends on the nature of 
the density field itself. We therefore normalise its values in the 
interval [0, 1] in order to apply it in a consistent way. 

There are other pos ible measures of feature inten- 



sity, 
trix 



(Fra ngi et al. , 1998b introduced the Frobenius 



ma- 



-^/Ij + A2 + A^ as a measure of second-order structure. 
However, this measure is biased towards blob-like structures 
and can produce erroneous signals in the detection of filaments 
and walls. 



6.5. Morphology Filter T 

For each level of the scale space, we can generate a 
Morphology Filter, T, from the Morphology Intensity Map I 
and Morphology Response Filter Ai. Formally we can write 
this as 

T = I<S>M (12) 

where the combination operator i8> simply means that every 
pixel of the Morphology Intensity Map, I, is multiplied by the 
value of the corresponding pixel in the Morphology Response 
Filter M. As described above, these hold information on dif- 
ferent aspects of the structural morphology, and by combin- 
ing them we can hope to improve on the results that would be 
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Fig. 8. The Scale Space Map Stack *P: the formalism selects at any location the scale with the optimal signal of the feature map. 
Depicted are the Feature maps T for three different scales (top row), and the resulting Map Stack ^ (bottom row), determined 
over the full range of scales. 



obtained by using either of them alone. Thus the Morphology 
Filter has its most significant values at those places where the 
morphology is close to what we are looking for. 



6.6. Feature Map f 

This is where, for each level of scale space, we combine infor- 
mation contained in the morphology mask S and filter T: we 
select out those regions of T where the morphology constraint 
is valid. 

For each level of the scale space, we can generate a Feature 
Map, 'F. The feature map comprises the information contained 
in the Morphology Filter 7~ and allocates it to the locations 
contained in the Morphology Mask S. Formally we can write 
this as 



(13) 



where the combination operator ® simply means that every 
pixel of the Morphology Filter, T, is modulated by the mask 
value £, 1 or dependent on whether the morphology con- 
straint is vaUd at the corresponding location. See figure 



6.7. Scale Space Map Stack ^ 

Each level of the scale space has its Feature Map constructed 
according to equation (fTSl ). They must now be combined in or- 
der to produce the definitive scale independent map of features. 
We can refer to *!* as the "feature stack" and formally write 
it as 



(14) 



levels n 



where the combination operator l+J represents a new Feature 
Map built by combining the individual Feature Maps, of 
the scale space. Each pixel of ^' takes on the maximum value 
of the corresponding pixel values in the stack of Feature Maps 
T^n in the Scale Space. We can write this (for a 3-D map) as 



^(i, j, k) — max 

Levels n 



r„{i,j,k) 



(15) 



where /, j, k represent the location of the pixels in the map. In 
this way we assign each point of the dataset a value quantifying 
the degree to which it can be said to be a part of some feature 
(blob, filament, or wall) on any of the scales investigated by the 
scale space. See figure (O. 
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(c) Threshold determination for walls 



6.8. Assigning Points to Features 

The Scale Space Map Stack ^' has to be thresholded in order to 
identify the most significant features. Tiis will be discussed in 
detail in Section O. It is at this point that we see astronomical 
input by requiring that the sought-after structure correspond to 
some structure that we would recognise. 

Given the Scale Space Map Stack for a given feature 
(blobs, filaments or walls), we can assign each particle of the 



Fig. 9. Thresholds for feature isolation based on the feature ero- 
sion criterion. The selected value is shown as a dotted vertical 
line. The object count to the right of the line declines due to 
erosion. 



original dataset to the specific feature identified in the Scale 
Space Map Stack. 



7. Cosmological Feature Detection: 
Threshold definition 

7.1. Texture Noise 

The final stage of each cycle of the analysis is the thresholding 
of the scale space map stack in order to identify individual ob- 
jects that are being sought in that cycle. Without the threshold- 
ing the maps are noisy and over-structured: we can refer to this 
as as "texture noise". This texture noise is simplest removed by 
applying a simple threshold to the processed maps. There is a 
potential problem in applying a simple threshold: it is neces- 
sary to determine a threshold that removes texture noise (how- 
ever that is determined) while leaving the sought-after features 
intact. 

7.2. Object Erosion Threshoid 

We set the thresholds for each feature to the value such that 
raising the threshold higher would start eroding objects and de- 
crease their number. In other words, the threshold value is set 
so that the object count is maximised while at the same time 
texture noise is eliminated. 

7.3. Identifying blobs 

We use tb to denote the value of above which a pixel is 
considered as part of a blob. Figure [9(a)] plots the number of 
objects detected above each value of the threshold, tb- 

For blob finding the thresholding is quite straightforward. 
At very low threshold, there will be many round objects (the 
eigenvalue criterion fixes their shape) of which only a small 
fraction will be the blobs we are seeking. As the threshold is 
raised from zero, the noise and the the less significant blobs are 
eliminated. There comes a point when the threshold stops an- 
nihilating these small, less significant, blobs and simply starts 
eroding the large blobs. This is the point where we define out 
optimal threshold. The dotted vertical line indicates the best 
value of Tb- 

If we plot a graph of the fraction of the sample volume oc- 
cupied by below-threshold blobs against the threshold we obvi- 
ously find a monotonic curve that rises from zero to one. This is 
shown in Figure [9(a)] where we see a two power-law behaviour 
with a break marking where the transition from texture noise 
annihilation to simple blob erosion takes place. 
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7.4. Identifying Filaments and Walls 

For filament and wall finding we again choose to threshold the 
distributions, but this time we decide on the optimal value of 
the threshold on the basis of the population curve of features 
defined at each threshold value. 

7.4.1 . Filaments 

We use Tf to denote the value of *P above which a pixel is con- 
sidered as part of a filament. Figure [9(b)| plots the normalised 
number of objects detected for each value of the threshold, Tf. 

The explanation for the shape of this curve is as follows. 
The low threshold (small-Tf ) objects are largely due to texture 
noise: the number of these declines as the threshold increases. 
When real filamentary features appear the number of detec- 
tions increases with Tf to reach a maximum. This is because 
at lower thresholds the features tend to percolate, so that rais- 
ing the threshold breaks the structure up into a greater number 
of filamentary objects. As the threshold rises further the fila- 
ments are eroded and get rarer. The point at which filament 
erosion starts to act is taken as the optimal value of Tf. This is 
indicated by the dotted line in the figure. 

7.4.2. Walls 

We use Tw to denote the value of *P above which a pixel is 
considered as part of a wall. Figure [9(c)| plots the normalised 
number of objects detected for each value of the threshold, t^. 

The threshold for defining walls is determined in the same 
way as for filaments. Note, however, that the particles classified 
as lying in blobs and filaments have been removed in previous 
cycles of the analysis so there is no longer a significant texture 
noise component. As the threshold is varied there is a peak in 
the number of walls that are found. At thresholds below this 
critical value the walls join up and percolate, eventually leaving 
one vast percolating structure. At higher threshold values walls 
are eroded and eventually destroyed. The dotted vertical line 
indicates the best value of tw 

7.5. Pseudo-code 

We have described the process of constructing a Feature Map 
and identifying features in that map. However there is a com- 
plication that arises in practise because both the Intensity Map 
and the Morphology Filter are built on a hierarchy of A values. 
In the case of the Morphology Filter, the different A's come 
in through equations (|9]l and (fTOb . In the case of the Intensity 
Map, different A's define the strength of different features as 
described in equation (fTTT i. 

The analysis cycle can be expressed in pseudo-code (see 
accompanying code in next column). In this form of pseudo- 
code, keywords (which correspond to class methods in object 
oriented programming) are in boldface. 

The nature of the hierarchy is such that we have first to 
identify blobs, remove them from the sample, then identify fil- 
aments, and after removing them from the sample finally iden- 
tify the walls. This arises because data points in blobs are de- 



fined by having three significant eigenvalues, data points in 
filaments are defined by having two significant eigenvalues, 
and data points in walls have only one significant eigenvalue. 
Identifying a filament before eliminating blobs would not work 
since the blobs would be more strongly detected. 



get PointSet 

set Feature = Blobs 

: Map_Feature 

resample PointSet to Mesh using DTFE 
construct ScaleSpace Hierarchy 

for each Level in ScaleSpace 

I 

build Hessian Eigenvalue Maps 

build using Eigenvalue Criteria for Feature 

1 

Morphology Mask, fi 
Feature Shape Fidelity, .S 
Morphology Response Filter, MiS) 
Feature Intensity Map, I 
I 

generate 

I 

Morphology Filter, T = I ® M 
Feature Map, T = 6 T 

1 

) 

stack ScaleSpace Feature Maps, 1* = l+J !F 

threshold Feature Maps using Feature Threshold Method 

in thresholded regions 

( 

identify Points 

publish Points 

remove Points from PointSet 

I 

if Feature = Blobs 

set Feature = Filaments 
else if Feature = Filaments 

set Feature = Walls 

else 

quit 

goto Map_Feature 



8. Areas of further development 

The methodology we have presented is very simple, yet, as we 
shall see, it is highly effective in differentiating the three main 
structural features that make up the cosmic web. The following 
section will test the methodology against a sample with con- 
trolled clustering: the Voronoi model, and present results for an 
N-Body simulation. Before going on to that analysis it is worth 
making a few remarks about some details of our procedure that 
might be enhanced. 

Our use of isotropic Gaussian filters is perhaps the most 
important limiting factor in this analysis. The largest filter ra- 
dius which is chosen is substantially smaller than the lengths 
of the typical filaments. Only the shorter filaments will get 
isotropised and they are "lost" since they make no contribu- 
tion in the scalespace stack. Our algorithm is indeed a long thin 
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filament finder. The main side-effect of the Gaussian smooth- 
ing is to make the profile (perpendicular to the filament) of the 
sharper (narrow) filaments Gaussian. A narrow filament having 
high density contrast will, under linear Gaussian smoothing, 
spill over into the large scales at a variety of thresholds and 
it will appear to be fatter than it really is. This latter problem 
is a consequence of applying simple linear filters: it is gener- 
ally overcome within the scal e space context by using nonlinear 



filters or by using wavelets jMartmez et al. 
I2007h 



2005 



Saar et al 



Another area for improvement is to use the eigenvectors as 
well as the eigenvalues themselves. Here we have simply relied 
on the relative magnitudes of the eigenvalues as indicators of 
curvature morphology. Had the eignevectors themselves been 
uncorrelated we might have concluded that there was structure 
when in fact there was only noise: the eigenvector coiTelations 
are good indicators of noise. 

A third area for improvement would be to use anisotropic 
smoothing filters. This leads us into another related approach to 
this problem: the use of nonlinear diffusion equations to locate 
structural features. This will be the subject of another article 
later on. 



9. Voronoi Clustering models 

To test and calibrate the Multiscale Morphology Filter we have 
applied the MMF to a set of four Voronoi Element Models. 
These models combine the spatial intricacies of the cosmic web 
with the virtues of a model that has a priori known proper- 
ties. They are particularly suited for studying systematic prop- 
erties of spatial galaxy distributions confined to one or more 
structural elements of nontrivial geometric spatial patterns. The 
Voronoi models offer flexible templates for cellular patterns, 
and they are easy to tune towards a particular spatial cellular 
morphology. In the case of the Voronoi models we have exact 
quantitative information on the location, geometry and identity 
of the spatial components against which we compare the out- 
come of the MMF analysis. 

9.1. Voronoi Models 

Voronoi Clustering Models are a class of heuris t ic models fo r 
cellular distributions of matter van de WeygaertI ( 1991 . 2002). 
They use the Voronoi tessellation as the skeleton of the cos- 
mic matter distribution, identifying the structural frame around 
which matter will gradually assem b le during the emer gence 
of cosmic structure IVoronoil (Il908r) : lOkabe et al. I (I2OOOI) . The 
interior of Voronoi cells correspond to voids and the Voronoi 
planes with sheets of galaxies. The edges delineating the rim 
of each wall are identified with the filaments in the galaxy dis- 
tribution. What is usually denoted as a flattened "supercluster" 
will comprise an assembly of various connecting walls in the 
Voronoi foam, as elongated "superclusters" of "filaments" will 
usually consist of a few coupled edges. The most outstanding 
structural elements are the vertices, coiTesponding to the very 
dense compact nodes within the cosmic web, rich clusters of 
galaxies. 



Model 


% blob 


% filament 


% wall 


% field 


A 


40 


30 


25 


5 


B 


43 


17 


32 


8 


C 


23 


37 


33 


7 


D 


27 


23 


42 


8 



Table 2. Voronoi Clustering Models. Percentage of galax- 
ies/points in the various morphological elements of the model. 



A more detailed description of the model construction may 
be found in Appendix IB. II We distinguish two different yet 
complementary approaches, Voronoi Element Models and kine- 
matic Voronoi models. 

Simple Voronoi models confine their galaxy distributions 
to one of the distinct structural components of a Voronoi tes- 
sellation: 

• Field 

Particles located in the interior of Voronoi cells (and thus 
randomly distributed across the entire model box) 

• Wall 

Particles within and around the Voronoi walls. 

• Filament 

Particles within and around the Voronoi edges. 

• Blobs 

Particles within and around the Voronoi vertices. 

Starting from a random initial distribution of points, these 
galaxies are projected onto the relevant wall, edge or vertex of 
the Voronoi cell in whose interior they are initially located. 

For our study we generated four different Voronoi cluster- 
ing models, labelled as A, B, C and D. They are all based upon 
a Voronoi tessellation generated by M = 53 nuclei distributed 
within a box of size L = 100 /i^'Mpc. The models are com- 
posite Voronoi Element Models and consist of the superposi- 
tion of galaxies located in field, walls, edges and vertices of 
a Voronoi tessellation. Our four test models contain - 32^ 
galaxies. The fraction of galaxies in the various components 
is a key parameter of the model, and is specified in Table 19.21 
In and around the walls, edges and vertices the galaxy distribu- 
tion follows a radial Gaussian density profile, with scale factors 
cTw = 1.0 /i^'Mpc, o-y- 1.0 /z^'Mpc and ctb - 0.5 /z^'Mpc. 

9.2. MMF Processing 

A considerable virtue of the Voronoi clustering models is that it 
is a priori known which galaxies reside in the various morpho- 
logical components of the Voronoi test models. This allows an 
evaluation of the absolute performance of the MMF and other 
morphology detection techniques by determining the fraction 
of the galaxies which are correctly identified as vertex, filament 
and wall galaxy. 

For each Voronoi model we computed the DTFE den- 
sity field from the particle distribution and applied the MMF. 
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Fig. 10. Recovered particles in Blobs, Filaments and Walls from a voronoi particle distribution. Particles inside blobs are detected 
(left), at 90/15 percent real/false detections. From the new blob-free distribution we detect particles in filaments (center) at 90/10 
percent real/false detections. Finally the blob-filament-free distribution is used to find the particles inside walls (right) at 80/10 
percent real/false detections." 



Following our previously described scheme, we first identified 
the blobs from the complete particle distribution. After removal 
of the blob particles, the filaments are found. Following the 
equivalent process for the filaments, the last step of the MMF 
procedure concerns the identification of the wall particles. The 
remaining particles are tagged as field particles. 



Figure ( fTOl i shows the outcome of the MMF applied to 
Voronoi Model C. Visually, the resemblance between real and 
MMF identified blob, filament and wall particles is remarkably 
good. The second row of panels shows the real detections of 
MMF: MMF clearly manages to identify all clusters, filaments 
and even the more tenuous walls in the weblike galaxy dis- 
tribution. The false detections do appear to have a somewhat 
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broader spatial distribution than those of the corresponding real 
detections. Most of them reside in the boundary regions of the 
blobs, filaments and walls: they are mainly an artefact due to 
the fact that the effective extent of the MMF morphology masks 
is slightly larger than the intrinsic extent of the Voronoi com- 
ponents. Fine-tuning of the filter scales (eq. [U is a potential 
solution for curing this artefact. 

9.3. Detection rate and Contamination 

The detection rate of blob, filament and walls galaxies is deter- 
mined and defined as follows. The galaxies in an MMF blob, 
filament or wall Map Stack ^' which are genuine Voronoi clus- 
ter, filament or wall galaxies are tagged as real detections. A 
galaxy detected by one of the three map stacks ^'t,, ^P/ or W^, 
intrinsically belonging to another morphological component is 
considered a false detection. For instance, a filament galaxy de- 
tected by '^h is a false blob galaxy. 

The main tunable parameters for optimizing the number of 
detected galaxies are blob, filament and wall threshold values, 
Tb, Tf and Th,. By lowering the blob threshold level ti,, defined 
through a regular density threshold (see sect. |7]), the number 
of MMF detected blob galaxies increases. The same holds for 
adjusting the filament and wall thresholds, in terms of the low- 
ering of the ^'f and levels. The galaxies detected by MMF 
include both real and false detections. As the threshold levels 
are adjusted the number of both will tend to increase. 

The detection rate at a given threshold level is the fraction 
of genuine blob, filament or wall galaxies which have been de- 
tected by the MMF. Ideally one would want to trace them all 
and have a 100% detection rate, in practice this is set by the the 
applied threshold. Based upon the 1-1 relation between Tb, t/ 
and T„ on the one hand and the corresponding blob, filament 
and wall detection rate on the other we use the detection rate as 
threshold parameter. 

The ratio of the corresponding number of false blob galax- 
ies to the total number of genuine blob galaxies is the blob 
contamination rate rate. The filament and wall contamination 
rate are defined in a similar way. Because a lowering of the 
threshold levels will result in a larger number of detections, 
both real and false, the contamination rate will be an increasing 
function of the detection rate. Note that the contamination rate 
may exceed 100% in the case the number of false detections 
exceeds that of the total number of genuine (blob, filament or 
wall) galaxies. 

9.4. Comparison 

We compare the MMF segmentation of the Voronoi models in 
blobs, filaments and walls with that achieved by a more di- 
rect criterion, that of a straightforward density threshold on the 
DTFE density field. We assign the label "DTC" to this naive 
procedure. 

Each of the morphological elements are identified with a 
particular (disjunct) range of density values. Blobs, ie. clusters, 
are identified with the highest density values. Filaments are as- 
sociated with more moderately high density values. Walls fol- 



low with density values hovering around unity to a few, while 
the field/voids may include densities down to a zero value. This 
approach has frequently been used to discriminate between 
virialized haloes and the surrounding matter distribution, and 
has even be en used in a n attem pt to define filamentary or pla- 
nar features lDolag et al. I (l2006h . However, it seriously oversim- 
plifies and distorts the perceived structure of the cosmic web. 
(This is presumably because filaments and walls differr in den- 
sity and have significant internal density structure. The simplis- 
tic density threshold approach does not reflect the reality of the 
structure: the range of densities in filaments overlaps with den- 
sities in walls and even with those of the outskirts of clusters. 



(Hahn et al.Ll2007h reach similar conclusions. 



9.5. Test results 

Figure [TT] compares the contamination rate as a function of 
the detection rate for the four different Voronoi models. The 
A,B, C and D models are distinguished by means of line style. 
The black lines relate to the MMF detections, the grey lines 
show the results of the equivalent DTC procedure. We find the 
following: 

• For all models, and for all morphologies, the MMF proce- 
dure is clearly superior to the DTC detections in suffering 
significantly lower contamination rates. 

• The MMF contamination is least for the blob detections. 
The filament contamination is lower than the wall contami- 
nation for models with many intrinsic filament galaxies (A 
and C). For models B and D, containing more wall galax- 
ies, the situation is the reverse. The same holds true for the 
DTC detections, be it much more pronounced and less fa- 
vorable wrt. the MMF detections. 

• The MMF and DTC blob contamination rate is more fa- 
vorable for the A and B models. Both models contain a 
relatively high fraction of blob galaxies. 

• The DTC blob contaminations are surprisingly bad, given 
that clusters are compact objects of high density with 
sharply defined boundaries. 

• The filament contamination rate is worse for models B and 
D, both marked by a relatively low amount of intrinsic fila- 
ment galaxies. This is true for both DTC and MMF. 

• The DTC contamination is extremely bad for models B and 
D, quickly exceeding 100%. This reflects the huge overlap 
in density range of filaments and other morphologies result- 
ing in a systematic inclusion of particles belonging to other 
morphologies. 

• For the MMF procedure there is a clear correlation be- 
tween the intrinsic wall galaxy fraction and the contamina- 
tion rate: model D has the highest number of wall galaxies 
and the lowest contamination. This is not true for DTC. 

In summary, we find that MMF clearly performs much 
better in tracing blob, filament and wall galaxies than a pure 
threshold criterion would allow. By comparing Voronoi mod- 
els A, B, C and D we find that MMF performs better for com- 
ponents which are relatively more prominent. Because of the 
mixture in densities between blobs, filaments and walls this is 
not necessarily true when using a simple density criterion. The 
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Fig. 11. Reals versus false detections for different voronoi models (see table |9.2| ) (A: solid, B:dotted, C:dashed, D:dotted-dashed) 
for blobs (left), filaments (center) and walls (right). We applied the MMF (black) and simple density thresholding (grey) in order 
to compare both methods. 



latter involves often excessive levels of contamination between 
galaxies in different morphological entities. If anything, this is 
perhaps the strongest argument for the use of the shape and 
morphology criteria enclosed in the MMF. 



10. N-body simulations 

The Large Scale Structure of the universe contains an intrincate 
mixture of morphologies. The boundaries separating each mor- 
phological component is rather ill-defined: clusters of galaxies 
are interconnected by filaments which in turn define the edges 
of walls. 

In order to explore the response of the MMF in this com- 
plex scenario we performed a cosmological N-body simulation. 
We give here only a few preliminary results to illustrate how the 
methodology works with a "real" galaxy distribution. A more 
detailed exploratin follows in a later paper. 

The simulation represents a LCDM model with Qa = 0.7, 
D.,,, = 0.3, h = 0.73 in a periodic box of side ISOMpc con- 
taining 256^ dark matter particles. We also run the same sim- 
ulation lowering the res olution to 128"^ part icles according to 
the prescription given by Kl vpin et al.l (120011) in order to assess 
the effect of mass resolution in the density field determination. 
For the scales studied here there is no significant difference be- 
tween the density fields computed from the two simulations 
since the mean interparticle separation is small enough to re- 
solve the intermedium-density structures ( Schaapl 2007). 



10.1. Results 

Figure [12] shows the result of applying the MMF to this sim- 
ulation. The multiscale nature of the MMF is clearly seen in 
figure[T2b which shows blobs of different sizes containing sim- 
ilar sized clusters of points (FigurjT2^). 



In the case of filaments and walls (see FigurdTSb andfTSlF) 
the multiscale nature of the MMF is less so obvious, however 
it is nonetheless there. 

It is clear from figure[T2]that even though the LSS presents 
a great challenge it can succesfully recover each morphological 
component at its characteristic scale. 



10.2. Blobs and Clusters 

Figure[T3K shows a section through the N-body model and fig- 
ure [T3b shows the blobs that are found in that section by the 
MMF process. Figure [T3t shows how these blobs relate to the 
underlying structure displayed in panel A. 

In an attempt to see how these blobs compare with clus- 
ters that would be identified via a more traditional method, we 
have use the HOP algorithm to find clusters. The HOP clus- 
ters are shown in figure [T3t) superposed on the MMF blobs. 
Superficially, we see that the agreement is remarkably good: 
the MMF blobs are indeed what we would subjectively call 
clusters. As in figure[T2]we can appreciate the scale range re- 
vealed by MMF. 

Making this comparison more precise is rather difficult ow- 
ing the the vastly different approaches to blob-finding. This will 
be discussed in detail in a subsequent paper dealing specifically 
with the application of MMF to N-body simulations. 

10.2.1. Filaments 

In Figure[T4]we show the particles that belong to the filaments 
defined at various scales of the scale space. The top left panel 
of figure [14] shows a histogram of the number of particles con- 
tained in the filaments seen at smoothing scales from 1 - 4/;"' 
Mpc. As expected the number of particles rises rapidly with 
smoothing scale (the filaments are fatter on larger scales and 
so encompass greater volume). The other three panels show 
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Fig. 12. MMF applied to N-body simulation. The top row shows a subsample (a) consisting of 10 % of the total number of parti- 
cles together with, in panels (b) and (c), the structures resulting from simple density thresholding using two different thresholds. 
Panels (b) and (c) both contain both spherical and elongated structures: there is a large amount of cross contamination between 
morphologies. Simple density thresholding is not an effective morphological discriminator. The second row shows the results of 
applying the MMF procedure showing clearly segregated (a) blobs, (b) filaments and (c) walls (for clarity we display only the 
largest structures. The thkd row shows the particles associated with the MMF defined structures. 



the points contained in filaments, seen in a slice of the N-body 
simulation at different resolutions. 

When these are stacked, application of equation ( fTSl l deter- 
mines whether a given pixel is a part of a filament. The process 
yields the filamentary map of figure [T2l 

10.3. Inventory of structures 

Finally we can simply count up the mass fraction of the model 
in various structural entities and the volume occupied by such 



structures to see how much of this N-body Universe lies in 
which structures. The result is shown in the pie diagrams of 
Figure [14] 

The result is hardly surprising: the clusters and filaments 
occupy about the same mass fraction and together contain more 
than half the haloes in the simulation. The clusters occupy by 
far least volume - they are dense systems and they are denser 
than the filaments. Recall, however, the important remark that 
we could not use density threshold alone to define these struc- 
tures (see the top row of panels in figure[T2]i. 
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Fig. 13. Comparing blobs found from HOP and from MMF. (A): Particles. (B): Isosurfaces of the blob identified with MMF. 
(C): Particles inside the blobs (black) and background particles in grey. (D): The position of the HOP Haloes (circles) and the 
particles inside the MMF blobs (dark grey). Light grey particles are just the rest. 



The large volume occupancy of filamentary structures ex- 1 1 . Conclusions and comments 

plains why our impression of the cosmic matter distribution is 
predominantly filamentary, and the fact that they are all long 
and thin (as illustrated in figure [T4l i emphasises the web-like 
nature of the structure. 



Perhaps the only surprise in this analysis is the relatively 
low volume occupancy of the walls in comparison with the fil- 
aments. This may be in part because most of the walls have 
broken up by the present epoch. It may also be in part due to the 
fact that the low number of particles in walls makes it relatively 
difficult to find them: they may get mis-classified as being part 
of the general field. It is difficult to assess this on the basis of 
the present experiments alone. 



MMF, our simple embodiment of Hessian-based scale space 
analysis of large scale cosmic structure, is remarkably sucess- 
ful in delineating the different structures that make up the cos- 
mic web. Since the morphology filters give us a direct mea- 
surement of blobness, filamentariness or wallness they can be 
used to characterize and quantify, in a systematic way, the large 
scale matter distribution. The technique has been tested using 
N-Body and Voronoi models. 
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Fig. 14. Particles defining filamentary structures in a slice of an N-body model. The grayscale images show the MMF detection 
of filamentary features on various filtering scales. Top lefthand: the filament volume occupancy (number of sample grid cells 
with a filament signal) as a function of smoothing scale. 



11.1. Void finding 

It should be emphaised that MMF is not a void finder except 
insofar as anything that is not in a blob, filament of wall might 
be deemed to be in a void region. In that case MMF would 
be a suitable tool for finding so-called "void galaxies" without 
being able to identify the host void. Void finding per se is al- 
most certainly best achieved via the Wat ershed (WVF) method 
Platen, van de Wevgaert & JonesI ( l2007b . 



11.2. Enhancements 

There are many areas where the MMF treatment could be en- 
hanced and some of these will be presented in future papers. 



We summarize a few issues here in order to place the present 
work in a more general perspective. 



The definition of the intensity component of the mor- 
phology filter could be improved by including other local 
propeties such as gradient, direction of eigenvectors, con- 
nectivity, etc. 

The Gaussian kernel is not the only possibility for produc- 
ing the scale-space representation: alternative kernels may 
improve the performance of the MMF. One side effect of 
using a simple Gaussian filter is that high peaks in high 
density filaments are detected always at larger scales even 
when their density profile is relatively narrow (the "filter 
smearing" we referred to earlier). 
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not be so easy! Bearing that in mind, the following is a list of 
possible application areas. 

The technique can readily be extended to analysis of veloc- 
ity data of various kinds such as Fingers Of God in cosmolog- 
ical redshift surveys, analysis of dynamical phase spaces, fea- 
ture detection in solar images, morphological characterization 
of structure in spiral arms, feature detection in radio datacubes, 
etc. Finding clusters and their substructures using MMF would 
provide an important alternative to HOP. Finding small, low 
surface brightness, galaxies in noisy neutral hydrogen surveys 
would be another useful application. 



(a) Volume occupied by each of the 
structural features. 

Mass content 




(b) Fraction of the mass occupancing 
each of the structural features. 



Fig. 15. Occupancy of Cosmic Web features, by volume (top) 
and mass(bottom) for a ACDM N-body simulation (see text). 

• Our implementation of the Multiscale Morphology Filters 
is grid-based and that required a resampling of the origi- 
nal point distribution data. It is possible to derive a similar 
set of filters using particle-based measures for the local dis- 
tribution of matter (e.g.: inertia tensor analysis), defining 
window functions and scale normalizations in a multiscale 
context. 

With respect to the above we would also like to refer to sect. [8] 
11.3. Applications 

The ability to accurately identify arbitrarily shaped struc- 
tures allows the possibility of seeking coiTelations within the 
structures that might otherwise be masked by other methods. 
Akeady, the method has been used to identify previoulsy un- 
known systemic p roperties in the alignment of haloes with the 



parent structures ( Arag on-Calvo et al.L 120071) . 



The technique has been illustrated in terms of spatial point 
data since that is relatively unambiguous. However, the MMF 
technique we have described is a quite general technique for 
scale-free feature finding: it only needs a mathematical pre- 
scription of what is being looked for, which in general may 
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Appendix A: 

The DTFE general reconstruction procedure 

For a detail ed specification of the DTFE density field procedure 
we refer to Schaapl ( 2007 !). In summary, the DTFE procedure 
for density field reconstruction from a discrete set of points 
consists of the following steps: 



• Point sample 

Given that the point sample is supposed to represent an un- 
biased reflection of the underlying density field, it needs to 
be a general Poisson process of the (supposed) underlying 
density field. 

• Boundary Conditions 

The boundary conditions will determine the Delaunay and 
Voronoi cells that overlap the boundary of the sample vol- 
ume. Dependent on the sample at hand, a variety of options 
exists: 

+ Empty boundary conditions: 

outside the sample volume there are no points. 
+ Periodic boundary conditions: 

the point sample is supposed to be repeated periodically 

in boundary boxes, defining a toroidal topology for the 

sample volume. 
+ Buffered boundary conditions: 

the sample volume box is suiTounded by a bufferzone 

filled with a synthetic point sample. 

• Delaunay Tessellation 

Construction of the Delaunay tessellation from the point 
sa mple. While we a l so stil l use the Voronoi-De launay code 
of lvan de WeygaertI (Il99lh and Ivan de Weygaert (,1994,) . at 
present there is a number of efficient library routines avail- 
able. Particularly noteworthy is the CGAL initiative, a large 
library of computational geometry routine^ 

• Field values point sample 

The estimate of the density at each sample point is the 
normalized inverse of the volume of its contiguous Voronoi 



- CGAL is a C++ library of algorithms and data structures for 
Computational Geometry, see www . cgal . org 
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cell 'Wj of each point /. The contiguous Voronoi cell of a 
point / is the union of all Delaunay tetrahedra of which 
point / forms one of the four vertices. We recognize two 
applicable situations: 

- uniform sampling process: 

the point sample is an unbiased sample of the underlying 
density field. Typical example is that of A^-body simulation 
particles. For D-dimensional space the density estimate is, 



p(x,) = (1+D) 



(A.l) 



with Wj the weight of sample point /, usually we assume 
the same "mass" for each point. 

- systematic non-uniform sampling process: 
sampling density according to specified selection process. 
The non-uniform sampling process is quantified by an a 
priori known selection function i^(x). This situation is typ- 
ical for galaxy surveys, i/'(x) may encapsulate differences 
in sampling density \l/{a,5) as function of sky position 
(a, S), as well as the radial redshift selection function i/'(r) 
for magnitude- or flux-limited surveys. For D-dimensional 
space the density estimate is , 



p(x,) = 



Wi 



(A.2) 



Field Gradient 

Calculation of the field gradient estimate V/|,„ in each 
D-dimensional Delaunay simplex m {D - 3: tetrahedron; 
D - 2: triangle) by solving the set of linear equations 
for the field values f at the positions r, of the (D + 1) 
tetrahedron vertices. 



For NN-interpol ati on see 
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Implementation of Natural neighbour 



interpolations is presently in progress. 
• Processing. 

Though basically of the same character, for practical pur- 
poses we make a distinction between straightforward pro- 
cessing steps concerning the production of images and sim- 
ple smoothing filtering operations and more complex post- 
processing. The latter are treated in the next item. Basic 
to the processing steps is the determination of field values 
following the interpolation procedure(s) outlined above. 
Straightforward "first line" field operations are Image re- 
construction and Smoothing/Filtering. 

+ Image reconstruction. 

For a set of image points, usually grid points, determine 
the image value, formally the average field value within 
the corresponding gridcell. In practice a few difl'erent 
strategies may be followed 

- Formal geometric approach 

- Monte Carlo approach 

- Singular interpolation approach 

The choice of strategy is mainly dictated by accuracy 
requirements. For WVF we use the Monte Carlo 
approach in which the grid density value is the average 
of the DTFE field values at a number of randomly 
sampled points within the grid cell. 

+ Smoothing and Filtering: 

A range of filtering operations is conceivable. Two of 
relevance to WVF are: 



/o /i fi h 



ro ri r2 rs 



(A.3) 



Evidently, linear interpolation for a field / is only meaning- 
ful when the field does not fluctuate strongly. 

Interpolation. 

The final basic step of the DTFE procedure is the field 
interpolation. The processing and postprocessing steps in- 
volve numerous interpolation calculations, for each of the 
involved locations x. Given a location x, the Delaunay tetra- 
hedron m in which it is embedded is determined. On the 
basis of the field gradient V/|„, the field value is computed 
by (linear) interpolation. 



/(X) = /(X,-) + V/|„,-(x-x,-). 



(A.4) 



In principle, higher-order interpolation procedures are also 
possible. Two relevant procedures are: 

Spline Interpolation 

Natural Neighbour Interpolation 



Linear filtering of the field / 
Convolution of the field / with a filter 
function Wj(x, y), usually user-specified. 



/.(X) 



/(x')W,(x',y)fl'x' 



(A.5) 



Natural Neighbour Rank-Ordered filterin g 
( Platen, van de Weygaert & JonesL 2007 ). 



• Post-processing. 

The real potential of DTFE fields may be found in sophis- 
ticated applications, tuned towards uncovering characteris- 
tics of the reconstructed fields. An important aspect of this 
involves the analysis of structures in the density field. The 
WVF formalism developed in this study is an obvious ex- 
ample. 

Appendix B: 

Voronoi clustering models 

Voronoi Clustering Models are a class of heuristic models fo r 
cellular distributions of matter dvan de Wevgaertlll99Ul2007h . 
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Fig. A.l. Three different patterns of Voronoi element galaxy distributions, shown in a 3-D cubic setting. The depicted spatial distributions 
correspond to a wall-dominated Voronoi Universe (top), a filamentary Voronoi Universe (centre) and a cluster-dominated Voronoi Universe 
(bottom). 



They use the Voronoi tessellation as the skeleton of the cos- 
mic matter distribution, identifying the structural frame around 
which matter will gradually assemble during the emergence 
of cosmic structure. The interior of Voronoi cells correspond 
to voids and the Voronoi planes with sheets of galaxies. The 
edges delineating the rim of each wall are identified with the 
filaments in the galaxy distribution. What is usually denoted as 
a flattened "supercluster" will comprise an assembly of various 
connecting walls in the Voronoi foam, as elongated "superclus- 
ters" of "filaments" will usually consist of a few coupled edges. 
The most outstanding structural elements are the vertices, cor- 
responding to the very dense compact nodes within the cosmic 
web, rich clusters of galaxies. 

We disti nguish tw o different yet complementary ap- 
proaches (see lvan de Wevgaert. 2007). One is the fully heuris- 
tic approach of "Voronoi Element models". They are particu- 
larly apt for studying systematic properties of spatial galaxy 
distributions confined to one or more structural elements of 
nontrivial geometric spatial patterns. The second, supplemen- 
tary, approach is that of the Voronoi Evolution models or 
Voronoi Kinematic models, which attempt to "simulate" we- 
blike galaxy distributions on the basis of simplified models of 
the evolution of the Megaparsec scale distribution. The Voronoi 
clustering models offer flexible templates for cellular patterns, 
and they are easy to tune towards a particular spatial cellu- 
lar morphology. To investigate the performance of MMF we 
use composite Voronoi Element Models, tailor-made heuristic 
"galaxy" distributions composed of a superposition of parti- 
cle distributions in and around the walls, edges and vertices of 
the Voronoi skeleton. A complete composite particle distribu- 
tion includes particles located in four distinct structural com- 
ponents: 

• Field 

Particles located in the interior of Voronoi cells (and thus 
randomly distributed across the entire model box) 

• Wall 

Particles within and around the Voronoi walls. 

• Filament 

Particles within and around the Voronoi edges. 



• Blobs 

Particles within and around the Voronoi vertices. 

For each of these four distinct distributions the model galaxies 
are projected onto the relevant Voronoi wall, Voronoi edge or 
Voronoi vertex or retained within the interior of the Voronoi 
cell in which they are located. Characteristic examples of sim- 
ple Voronoi Element galaxy distributions are the ones shown 
in the boxes in Fig. lA.ll The depicted distributions concern 
a wall-dominated Voronoi Universe (lefthand), a filamentary 
Voronoi Universe (centre) and a cluster-dominated Voronoi 
Universe (righthand). 

In the case of composite models the fraction of field galax- 
ies Xc, wall galaxies X„, filaments galaxies Xf and blob galax- 
ies Xb, with the constraint Xc + X^ + Xf -\- X/, - 100, is a key 
input parameter of the model. 

B.1. Initial Conditions 

The initial conditions for the Voronoi galaxy distribution are: 

• Distribution of M nuclei, expansion centres, within the sim- 
ulation volume V. The location of nucleus m is y,„. 

• Generate model galaxies whose initial locations, x„o (« - 
1 , . . . , A^), are randomly distributed throughout the sample 
volume V. 

• Of each model galaxy n determine the Voronoi cell 'Va in 
which it is located, ie. determine the closest nucleus ja. 

All different Voronoi models are based upon the displacement 
of a sample of A^ "model galaxies". The initial spatial distribu- 
tion of these A^ galaxies within the sample volume V is purely 
random, their initial locations x„o (n - I, . . .,N) defined by a 
homogeneous Poisson process. A set of M nuclei within the 
volume V corresponds to the cell centres, or expansion centres 
driving the evolving matter distribution. The nuclei have loca- 
tions y,„ (m = 1 , . . . , M). 

Following the specification of the initial positions of all 
galaxies, the second stage of the procedure consists of the cal- 
culation of the complete Voronoi track for each galaxy n - 
1, . . .,N (sec. IB. 21 ) towards its wall, filament or vertex, or its 
location within a cell when it is a field galaxy. . 
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Neighbouring 
nucleus * 



Fig.B.l. Schematic illustration of the galaxy projections in the 
Voronoi clustering model. See text. 



Simple Voronoi Element Models place all model galaxies in 
either walls, edges or vertices. The versatility of the model also 
allows combinations of element models, in which field (cell), 
wall, filament and vertex distributions are superimposed. The 
characteristics of the patterns and spatial distribution in these 
Mixed Voronoi Element Models can be varied and tuned accord- 
ing to the fractions of wall galaxies, filament galaxies, vertex 
and field galaxies. 

B.2. Voronoi Tracks 

The first step of the formalism is the determination for each 
galaxy n the Voronoi cell 'Va in which it is initially located, ie. 
finding the nucleus ja which is closest to the galaxies' initial 
position x„o- 

In the second step the galaxy n is moved from its initial 
position x„o, towards its final destination in a wall, filament or 
vertex (see fig. IB. 11 1. The first section of the galaxy displace- 
ment is the radial path along the direction defined by the galax- 
ies' initial location wrt. its expansion centre ja- This direction 
is defined by the unity vector e„a. 

If the galaxy is a field galaxy it remains at its original lo- 
cation. If it is a wall galaxy it is projected along direction e„a 
onto the Voronoi wall with which the radial path first in- 
tersects. Filament galaxies are moved along the wall to the lo- 
cation where the path intersects Aapy. Finally, if it is cluster 
galaxy the galaxies' path is continued along the edge Kapy until 
it reaches its final destination, vertex Safiys. The identity of the 
neighbouring nuclei ja, jf}, jy and jg, and therefore the identity 
of the cell 'Va, the wall S^,^, the edge Aa^y and the vertex H^^j^^, 
depends on the initial location x„o of the galaxy, the position y^, 
of its closest nucleus and the definition of the galaxies' path 
within the Voronoi skeleton. 



In summary, the path x„ is codified by 

(B.l) 

in which the four different components follow the directions 
defined by: 

unity vector of path within Voronoi cell 'Va 
unity vector of path within Voronoi wall Y.ap 

* ^nafiy - 

unity vector of path along Voronoi edge Aa^y 

• Vertex Ea/}ys 

The cosmic matter distribution is obtained by calculating 
the individual displacement factors (s„„, s„a/}, s„apy) for each 
model galaxy, corresponding to their location within either 
wall, filament or vertex. In the Voronoi Element models all 
galaxies are directly projected onto wall, edge or vertex follow- 
ing the path depicted in fig. IB. II The coiTesponding displace- 
ment factors in eqn. IB.2l for a wall, filament or cluster galaxy 
are 

Walls {s„a,S„al},S„apy) = (l'„,0,0) 

Filaments {s„a, Snap, s„apy) = {v„,o-n,0) (B.2) 

Clusters (s„a. Snap, SnaPy) — {l'n,Crn,A„) 

where the values of the parameters v„, o-„ and A„ characterize 
the crossing of the galaxies' path with the wall, edge or vertex 
towards which it moves. 

A finite thickness is assigned to all Voronoi structural ele- 
ments. The walls, filaments and vertices are assumed to have a 
Gaussian radial density distribution specified by the widths 
of the walls, Rp of the filaments and Ry of the vertices. Voronoi 
wall galaxies are displaced according to the specified Gaussian 
density profile in the direction perpendicular to their wall. A 
similar procedure is followed for the Voronoi filament galaxies 
and the Voronoi vertex galaxies. As a result the vertices stand 
out as three-dimensional Gaussian peaks. 



